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Abstract 

We present a new method for component separation from multifrequency maps 
aimed to extract Sunyaev-Zeldovich (SZ) galaxy clusters from Cosmic Microwave 
Background (CMB) experiments. This method is best suited to recover non-Gaussian, 
spatially localized and sparse signals. We apply our method on simulated maps of 
the ACBAR experiment. We find that this method improves the reconstruction of 
the integrated y parameter by a factor of three with respect to the Wiener filter 
case. Moreover, the scatter associated with the reconstruction is reduced by 30 per 
cent. 
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1 Introduction 


The study of the Cosmic Microwave Background (CMB) has greatly improved 
our understanding of the Universe in the last decade. The measurement and 
interpretation of the CMB power spectrum has allowed to determine the 
most important cosmological parameters with very high accuracy. More ex¬ 
periments, now planned or undeway, will produce higher resolution multi¬ 
frequency maps of the sky in the 100-400 GHz frequency range. One of the 
most important new scientihc goals of these experiments is the detection of 
clusters through their characteristic Sunyaev-Zeldovich (SZ) signature (Sun- 
yaev & Zeldovich 1980). Because the SZ signal is substantially independent 
of redshift, SZ clusters will be observed at very large distances, quite inde¬ 
pendently from their mass. Such clusters may be used to infer cosmologi¬ 
cal information via number counts and power spectrum analysis of SZ maps. 
Many studies have shown the great potential of these new technique. These 
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Table 1 

The characteristics of the ACBAR experiment. 

V (GHz) FWHM (arcmin) noise (;uK/beam) 


150 

5 

6 

220 

5 

10 

280 

5 

12 


estimates, however, typically assume that all clusters above a certain flux are 
perfectly reconstructed and detected in the CMB maps. In practice, this may 
not be the case. SZ clusters have radio intensities comparable to other in¬ 
tervening cosmological signals like the CMB and point sources. Despite the 
different frequency and spatial dependence of these signals, it is not so easy 
to disentangle them. Moreover, beam smearing and instrumental noise play a 
role in our ability to adequately reconstruct the observed cluster. These ar¬ 
guments raise the necessity to assess how well a certain technique performs 
in reconstructing the cluster signal given the experiment specihcations. For a 
given cluster Compton parameter y, the reconstructed value may also depend 
on the cluster location and shape. Therefore, there is an error associated with 
the reconstruction technique which needs to be assessed and accounted for 
when it comes to relate cluster’s observables with cosmological models. More¬ 
over, the specific observable to use may depend on the type of experiment 
in hand. In this paper we address these issues for the ACBAR experiment 
(see Pierpaoli et al. 2004 for a similar analysis applied to Planck and ACT). 
Several techniques have been developed for image reconstruction in the multi- 
component case (Herranz et al. 2002, Stolyarov et al. 2002). In most cases, 
these techniques are optimal in reconstructing the CMB fluctuation signal, 
which is Gaussian and well characterized in Fourier space. Clusters of galaxies 
maps present very different features from CMB ones, in particular: %) clusters 
are “rare” objects in the map, they don’t £11 out the majority of the space; 
ii) The cluster signal is non-Gaussian on several scales, and in particular on 
scales associated with the typical core size; in) Different scales are correlated 
in Fourier space. Keeping these characteristics in mind, in this paper we de¬ 
velop a method aimed to better reconstruct the SZ galaxy cluster signal from 
multifrequency maps. Our map reconstruction method is wavelet based and is 
best suited to reconstruct the specific non-Gaussian signal expected in galaxy 
clusters maps. In this work we use the simulated cluster maps by Martin White 
available at: http://pacl.berkeley.edu/tSZ/. These maps are not full hydrody- 
namical simulations, the gas here has been introduced using the dark matter 
as a tracer. However for the kind of experiment in hand which would not re¬ 
solve the cluster structure anyway, this should not present a major problem. 
The underlying cosmology corresponds to the concordance model. We used 
ten maps of 10 x 10 degrees. 
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2 The reconstruction method 


2.1 The wavelet decomposition used 

We use a two-dimensional overcomplete wavelet representation which is more 
adequate to the analysis of astrophysical images. The wavelet decomposition 
of a signal s in our case reads: 

J M 

S H II (1) 

q£N^ j=0 m=l qe2-J 

where (pq are the scaling functions, 'pj,m,q are the wavelets and (,) are scalar 
products. The sum 4'q)4>q is the projection of s on the coarsest scale, i.e. 
a low-pass version of s. Each scaling coefficient (s, pq) contains information 
about the signal s at the coarsest scale and at a specific location in space q. For 
j fixed, the sum 'Pj,m,q)'Pj,m,q IS the projection of s on the scale j, i.e. a 

band-pass version of s. Each wavelet coefficient (s, 'pj,m,q) contains information 
about the signal s at the specific scale j, orientation m and location in space 
q. As usual with wavelet transforms, changing scale is done by dilating, and 
changing location is done by translating the wavelet: pj^i^rn^qix) = 
and 'pj^rn,q{x) = 'pj^rnfl{x — q) Hence, scale j + 1 corresponds to a frequency band 
that is twice as wide and for which the central frequencies are twice as large as 
that of scale j. On the other hand, in space, the wavelets at scale j+l are better 
localized than at scale j since they are more narrowly concentrated around 
their center q (see Fig.l, column 1 and 2). Unlike the 2-D Daubechies wavelets, 
the wavelets (and scaling function) we use here are defined in the Fourier 
plane. This ensures that they are well concentrated in frequency. Moreover, it 
enables us to introduce orientation by rotating the Fourier transform of the 
wavelet (see Fig.l, column 2 and 3). If / is the Fourier transform of / and 
(r, 6) are polar coordinates, then: ‘pj^m^qi.f, 0) = pjp^qir, 0 — ^) The transform 
is therefore close to rotation invariant and computation is fast via FFT. The 
Fourier transform of the wavelets and scaling functions read: 


TT 

L(r) = cos (- log 2 (r))hi <,.<2 -h 5r<i 

low pass 

(2) 

TT 

H{r) = sin (- log 2 (r))(ji <^<2 -h Sr >2 

high pass 

(3) 

n iQ\ (M —1)! q\M-i 

Gm{,0)= ^- — |2cos0| 

^M[2{M -l)]\ 

oriented 

(4) 

MT,e) = L{2T,d) 


(5) 
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i>o, 
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(6) 
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Fig. 1. Top row: wavelets in space; Bottom row; wavelets in Fourier plane. Left to 
right; wavelet at a fine scale j + 1, centered at location qo, oriented along the first 
diagonal; wavelet at a coarser scale j, centered at location qi, oriented along the 
first diagonal; wavelet at the same coarser scale j, centered at location q 2 , oriented 
along the horizontal axis; scaling function, centered at location q 2 . 

The set of all wavelets and scaling functions determines a redundant system 
(they are linearly dependent), however, the Plancherel equation holds. 


2.2 The estimator 


Formally, our goal is to estimate several processes (CMB, SZ, point sources) 
from their contributions in observations at different frequencies. We estimate 
the processes {x{p,uq)}p from the observations {y{r')}u given that: y{v) = 
J2p f{P) ^)x{p, i^o) * B(iy) + where x{p, uq) is the template of the pth pro¬ 
cess at a given frequency uq, f{p, u) is the frequency dependence of the pth pro¬ 
cess, B{v) is the beam and N{v) is the frequency dependent white noise. Our 
estimation method will rely on two principles to discriminate the contributions 
from different processes. The hrst one is that we know some statistical prop¬ 
erties of the processes (e.g. the CMB and noise are Gaussian processes, while 
the clusters are not). The second one is that some spatial properties of the 
processes can be captured by modeling the coherence of wavelet coefficients. 
For example, clusters can be described as spatially localized structures with 
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a high intensity peak. To estimate a particular wavelet coefficient one 

describes the statistics of a neighborhood of coefficients around it by a Gaus¬ 
sian scale mixture. For example, = {Xj^rn,q,Xj^rn,q+l,Xj^rn,q-l,Xj_i^m,q) is 

a neighborhood of coefficients around Xj^m,q- If contains wavelet coefficients at 
the same scale with close location, and at a close scale with the same location. 
The Gaussian scale mixture is the model: 

X = ^/zu (7) 


where u is a centered Gaussian vector of the same covariance as x, the mul¬ 
tiplier 2 : is a scalar random variable and the equality holds in distribution, u 
and 2: are independent and E{z} = 1. The covariance of x captures the spatial 
coherence of the process. The (non-)Gaussianity of the signal is captured by 
the distribution of the multiplier 2:. To illustrate the idea for the reconstruction 
process, let us consider the simple case where we observe one process polluted 
by noise: y = x + N {in wavelet space), x is a Gaussian mixture x = \/z\i, 
with u a Gaussian vector. If 2; was a constant, z = Zq, then T’jxly, 2; = Zq), 
the Bayes Least square estimate of '^j,m,q given the observed vector yj,m,q and 
2;, would be the Wiener filter: 

F;{x|y, 2; = 2:0} = 2 :oCx(^;oCx + Cn)"V ( 8 ) 


where Cv is the covariance matrix of the vector v. However in our model, 
is not a constant, so F’jxly}, the Bayes Least square estimate of is a 

weighted average of the Wiener filters above: 

00 

^{x|y} = jp{z = 2 :o|y)£^{x|y, ^ = zo} dzQ (9) 

0 


The weights are determined by the probability of 2 ; given the observation 
yj,m,q, noted p{z = zo\yj,m,q), which is computed via Bayes rule: 


P{z = 2:o|y) 


P{y\z = Zo)p,{zo) 
fp{y\z = z')p^{z')dz' 


( 10 ) 


where p(y| 2 ; = z') is a centered Gaussian vector of covariance 2;'Cx -l- Cn, 
and pz is the probability distribution of 2 ;. Following this procedure, one gets 
an estimate i7{xj „j,g|yj,m,g} for each neighborhood of coefficients From 

this estimated vector, we only keep the estimate of central coefficient Xj^rn,q- 

In the case where p{z) = Sd{z — 0), the Gaussian scale mixture described 
in eq. (7) reduces to a Gaussian process, which is an accurate model for the 
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CMB signal. Other signals, in particular the cluster signal, are typically non- 
Gaussian. In order to model them, we will need a more elaborate distribution 
for In this paper we use a distribution p{z) that we derived from the input 
SZ maps with the technique described in Pierpaoli et ah (2004). The cluster’s 
distribution p{z) has a tail for high values which is caused by the high 
intensity points in the cluster centers. By using this distribution instead of 
the delta function (which would correspond to a Gaussian process) we are 
suggesting to the reconstruction method that in the map there should be 
more “high intensity” points than in the corresponding Gaussian case with 
the same variance. In Pierpaoli et ah (2004) we also describe other choices 
for p{z) and conclude that the hnal performance in reconstructing the cluster 
center doesn’t depend on the specihc shape of the distribution p{z) provided 
that p{z) has enough power in the “high-intensity” tail. 


3 Results 


We applied the method described above to AGBAR simulated maps with 
realistic beam and noise properties (see table 1). We ignored here the com¬ 
plications in the noise correlation matrix, and assumed instead white noise. 
We don’t expect this approximation to highly compromise the overall per¬ 
formance. In hgures 2 and 3 we show the performances of our method in 
reconstructing the y parameter of the largest and most intense clusters in the 
simulations. The reason for such distinction is the following: the method pro¬ 
posed here make use of the spacial covariances of the cluster’s signal. Therefore 
if the beam size of the experiment is comparable to the typical cluster size (as 
is the case with AGBAR) clusters with similar intensities but different core 
size may not be reconstructed equally well, since compact clusters are more 
likely to be confused with noise. The actual performance in the reconstruction 
largely depends on the noise level. For this reason we perform the two analysis 
for “extended” and “intense” clusters separately. 

We adopt the following procedure: hrst we compute the input and output y 
parameter integrated over a given angle (specihed on the x axis); we then 
fit a line to the input/output values and calculate the average departure of 
the output values from that line. The figures report the slope of the htted 
line (solid lines) and the scatter about it (dashed lines) for the p{z) cluster 
distribution (blue) and for a delta function corresponding to an hypothetical 
Gaussian signal (pink). 

In these hgures we notice that the Wiener hltering (which assumes the delta 
function for the p{z) distribution) underestimates the high intensity peaks 
corresponding to massive clusters. The estimator which uses the true p{z) 
distribution performs on average a factor of three better in reconstructing 
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ACBAR: y^g^>5e-05, Vj g>4.7e-05, yg>4.3e-05 



Fig. 2. The integrated output y parameter versus the input one for the (12) largest 
clusters in the simulations. These values are obtained by using the ACBAR specih- 
cations. The title indicates the y values of the less intense cluster considered: ymax, 
1 / 3.6 and 1/6 are the central intensities in the input map with no smoothing (the pixel 
size of the simulation being 1.17 arcmin) and with a smootihng angle of 3.6 and 6 
arcmin respectively. 


ACBAR: y|^gj|>2.9e-05, y3g>2.7e-05, yg>2.2e-05 



Fig. 3. The integrated output y parameter versus the input one for the (34) most 
intense clusters in the reconstructed maps. The values at the top have the same 
meaning as in fig. 2. These values are obtained by using the ACBAR specifications. 

the central intensity of the clusters, quite independently from the integration 
angle assumd. The error average departure of the output value from the htted 
line is also reduced by about 30 per cent. We conclude that the inclusion of 
non-Gaussian information greatly improves the reconstruction of clusters SZ 
signal from observed maps. Performances, however, depends on the specihc 
characteristics of the experiment in hand, as well as on the characteristics of 
the clusters that we aim to reconstruct (i.e. total intensity and dimensions). In 
the specihc ACBAR case, by comparing hgs 2 and 3 we notice that extended 
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clusters are on average slightly better reconstructed and present a significantly 
lower scatter. 
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